Discrete diffraction and shape-invariant beams in optical waveguide arrays 
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General properties of linear propagation of discretized light in homogeneous and curved waveguide 
arrays are comprehensively investigated and compared to those of paraxial diffraction in continuous 
media. In particular, general laws describing beam spreading, beam decay and discrete far-field 
patterns in homogeneous arrays are derived using the method of moments and the steepest descend 
method. In curved arrays, the method of moments is extended to describe evolution of global beam 
parameters. A family of beams which propagate in curved arrays maintaining their functional shape 
-referred to as discrete Bessel beams- is also introduced. Propagation of discrete Bessel beams in 
waveguide arrays is simply described by the evolution of a complex q parameter similar to the 
complex q parameter used for Gaussian beams in continuous lensguide media. A few applications of 
the q parameter formalism are discussed, including beam collimation and polygonal optical Bloch 
oscillations. 
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I. INTRODUCTION 



Linear and nonlinear propagation of 'discretized' light 
in arrays of evanescently-coupled optical waveguides has 
received a great and increasing interest in the past recent 
years (see, for instance, 0, 0] and references therein). 
As compared to diffraction or refraction in continuous 
(non-structured) media, discrete diffraction and refrac- 
tion in waveguide arrays show rather uncommon effects 
which result from the evanescent coupling among ad- 
jacent waveguides forming a one-dimensional or a two- 
dimensional lattice. For instance, linear propagation of 
light waves in homogeneous arrays may show diffrac- 
tion reversal and self-collimation effects [3,Jj| , anomalous 
refraction the discrete Talbot effect @, and quasi- 
incoherent propagation [6j to name a few. Remarkably, 
discrete diffraction can be tailored by properly introduc- 
ing inhomogeneities in the lattice or by varying its topol- 
ogy. In particular, since the first proposals and demon- 
strations of optical Bloch oscillations 0-Q and 'diffrac- 
tion management' in zig-zag arrays the use of waveg- 
uide arrays with curved optical axis has been extensively 
investigated both theoretically and experimentally, with 
the demonstration of diffraction suppression via Bloch 
oscillations [l|, [7rll0| or dynamic localization [l2|, fl3j . 
polychromatic diffraction management |14j . astigmatic 
diffraction control [Hj], multicolor Talbot effect [lj|, and 
discrete soliton management [l6j . Linear and nonlinear 
light propagation at the surface or at the interface of 
two waveguide lattices also exhibits a variety of inter- 
esting properties which have been investigated in several 
recent works (see, for instance, @, [l7l - [i"9j and references 
therein). In spite of such a great amount of works, some 
facets of discrete diffraction, even in the simplest linear 
propagation regime, have been overlooked. Though in 
the linear regime the impulse response (Green function) 
of the array may be rather generally calculated analyti- 
cally -either in straight or curved geometries and in pres- 
ence or not of boundaries- and its knowledge is enough 



to predict light evolution for any assig ned initial exci- 
tation condition (see, for instance, [l3l Il7j). some gen- 
eral issues of discrete diffraction, which are well known 
for paraxial propagation of beams in continuous media, 
have not been comprehensively addressed, including: (i) 
a description of global beam parameter evolution in a 
closed analytical form; (ii) far-field discrete diffraction in 
homogeneous array (the analogue of Fraunhofer diffrac- 
tion in homogeneous continuous media); (iii) the gen- 
eral scaling law of beam broadening and beam decay, 
especially close to the self-collimation condition (also re- 
ferred to as sub-diffraction) which is commonplace to the 
more general class of photonic crystal structures (see, 
for instance, [13]); (iv) the existence of shape- invariant 
discretized beams, i.e. special families of field distribu- 
tions which -like Gaussian beams in continuous lensguide 
media- do propagate in straight or curved waveguide ar- 
rays maintaining their functional shape. It is the aim of 
this work to shed some light into such issues. In par- 
ticular, it is shown rather generally that: (i) the scaling 
law describing broadening of discretized light in homo- 
geneous arrays is the same as that of standard parax- 
ial diffraction theory of homogeneous continuous media 
(beam size asymptotically grows linearly with propaga- 
tion distance), independently of the precise array disper- 
sion curve and even along self-collimation directions; (ii) 
in a homogenous array, the discrete far-field pattern is 
not the (discrete) Fourier transform of the near-field dis- 
tribution, and the scaling law of beam decay may depend 
on the observation angle; (iii) special field distributions, 
which propagate in straight or curved waveguide arrays 
maintaining their functional shape and referred to as 'dis- 
crete Bessel beams', can be introduced for simple tight- 
binding waveguide models; (iv) a discrete Bessel beam is 
defined by a complex q parameter, analogous to the one 
used for Gaussian beams in continuous lensguide media, 
and propagation of the q parameter along the array ad- 
mits of a simple geometric interpretation. 
The paper is organized as follows. In Sec. II general prop- 
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erties of discrete diffraction in homogeneous waveguide 
arrays are presented, including the derivation of the gen- 
eral scaling laws of beam broadening and beam decay, far- 
field discrete diffraction, with a a note on self-collimation 
regimes. In Sec. Ill, some general rules of beam prop- 
agation in curved waveguide arrays are derived within 
the nearest-neighbor coupling approximation, whereas 
in Sec. IV the family of shape-invariant discrete Bessel 
beams is introduced, together with the complex q param- 
eter formalism. Applications to beam collimation and 
polygonal optical Bloch oscillations are also presented. 
Finally, in Sec.V the main conclusions are outlined. 

II. DISCRETE DIFFRACTION IN A 
HOMOGENEOUS WAVEGUIDE ARRAY 

A. Continuous model of discrete diffraction 

The starting point of our analysis is provided by a 
rather standard model describing linear propagation of 
monochromatic light waves along the z direction of a 
one-dimensional or two-dimensional array of waveguides 
in the single band and tight-binding approximations. For 
instance, in a one-dimensional array such conditions are 
satisfied when the tilt of beams and waveguides at the in- 
put facet is less than the Bragg angle, so that the lowest- 
order band of the array is excited and beam propagation 
is primarily characterized by coupling between the funda- 
mental modes of the waveguides. For a two-dimensional 
array, the relevant equations describing discrete diffrac- 
tion in a single band approximation read 

iCn,m ^ ^7i-I.m-rQ,r (1) 

l,r 

where c„ iTO (z) is the complex amplitude of the fundamen- 
tal waveguide mode at the lattice site v n ^ m = na + mh 
identified by the indices (n,m), a and b are the lattice 
vectors of the unit cell, the dot denotes the derivative 
with respect to z, and A„ jm = AJ^ n are the coupling 
rates. In order to derive a general rule of beam broaden- 
ing due to discrete diffraction, it is worth introducing a 
continuous field envelope ip(x,y,z) satisfying the scalar 
Schrodingcr-likc equation 

id z iP(r,z) = H (p)iP(r,z), (2) 
where r = (x,y), p = — iV T , 

H (p) = A„ !tn exp (-ir„ !m • p) , (3) 

and r„ jm = na + mh. Taking into account that 
exp(— iR ■ p)tp(r,z) = f/i(r + R,z), it follows that the 
solution Cn >m {z) to the discrete equation (1) can be iden- 
tified with ip(r = na. + mh, z). The formulation of the 
discrete light propagation problem [Eq.(l)] as a contin- 
uous problem [Eq.(3)] is a well-established procedure in 



solid-state physics [2l| which enable the use of certain 
analytical techniques, such as the method of moments, 
developed for the continuous Schrodinger equation or for 
the paraxial wave equation (see, for instance, [22|,[23j]). In 
addition, the continuous model includes, as a particular 
case, the problem of paraxial diffraction in a homoge- 
neous medium (e.g. in the vacuum), which is attained 
by simply assuming for the Hamiltonian Hq(p), in place 
of Eq.(3), the (normalized) parabolic form 

H (P) = y- (4) 

The normalization conditions J dr\tp(r, z)\ 2 = I for 
Eq.(2), and Y, n , m \ c n, m (z)\ 2 = 1 for the discrete prob- 
lem (1), will be assumed in the following analysis. 

B. General law for beam spreading: moment 
analysis 

Two global parameters describing beam propagation 
are the beam center of mass (r) = (x)u x + (y)u y and the 
transverse beam spot sizes w x (z) and w y (z) defined by 

(r) = / dri#(r,z)| 2 , (5) 



w x (z) = ^/(x 2 ) - (x) 2 (6) 

w y (z) = vV> - (y) 2 (7) 

where 

(x 2 )(z) = J drx 2 \^r,z)\ 2 , (y 2 )(z) = J dvy 2 \^{v , z)\ 2 . 

(8) 

Note that the above definitions hold for both continu- 
ous and discrete diffraction models. In the latter case, 
assuming ?/>(r, z) to be a piecewise constant function in 
each cell of the lattice and taking |a x b| = 1 for the area 
of the unit cell, the integral over r may be replaced by 
a double sum over the cell indices n and m, i.e. in the 
discrete model one has / dr —> J2 m n . 
The evolution equations for r and w XtV can be readily 
obtained in a closed form by writing a set of Eherenfest 
equations for the expectation values of r, x 2 , y 2 , and of 
commutator operators that arise in the calculation. The 
expectation value (A) = J drip* (r, z)A(r, p)ip(r, z) of any 
operator A(r, p) (not necessarily self-adjoint) evolves ac- 
cording to 

^ = -i([A,H }) (9) 
and the following commutation relations 

[r,/(p)]=iV p /, [ff(r),p]=*V rS (10) 
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hold for any functions /(p) and g(r). For A = r, one 
obtains 



d(r) 

dz 



= (Vp/fo), 



d(V p H 0/ 



dz 



= 0, 



i.e. 



(r)(,)Hr)(0) + (V p % 



(11) 



(12) 



which is the evolution equation of the beam center of 
mass. Note that the path followed by any beam is al- 
ways straight, regardless of the specific form of H or 
initial field distribution which just determine the trans- 
verse drift velocity (V p 7Jo) of the beam. To determine 
the evolution equation of the beam spot size w x , let us as- 
sume A = x 2 , so that the following cascade of Eherenfest 
equations [Eq.(9)] is obtained 



d dH 

dz dp x 
d ( 
dz 



d(x 2 ) 
dz 
dH , 



{x dHo + dH 
dp x dp x 

2 



dp, 
dp x 



-x 



dp x 



(13) 
(14) 
(15) 

where the expectation values on the right hand side of 
Eq.(16) are calculated at z = 0, i.e. for the initial beam 
distribution. From Eqs.(6), (12) and (16) the following 
evolution equation for the beam spot size w x is then ob- 
tained 



After integration, one obtains 

{ x 2 )(z) = (x 2 )(0) + (x^ + ^x)z + ( 

OPx OPx 



dH 

dp x 



w 2 (z) = w 2 (0) + a x z + (3 2 z 2 , 
where we have set 

a x = ((x - (x))^^ + *^-{x - (x))), 



dp x dp x 



,{dH 



V dp* 



,ggo. 
' dp x ' 



(17) 



(18) 



(19) 



and the expectation values are calculated at z = 0. A 
similar expression for w y (z) can be obtained by replac- 
ing x with y in Eqs.(17), (18) and (19). A major re- 
sult expressed by Eq.(17) is that, regardless of the par- 
ticular form of £/o j w x (z) (and similarly w y (z)) asymp- 
totically grows with z linearly, with a diffraction length 
given by ~ l/Ar- Therefore -and this one of the ma- 
jor result of this section- the broadening law of a spatial 
beam due to diffraction does not differ for discrete or 
continuous diffraction. In addition, for a beam carry- 
ing a finite power and admitting of finite moments (a; 2 ) 
and (y 2 ), the coefficient (3 2 given by Eq.(19) is always 



strictly positive and does not vanish. This can be gen- 
erally proven by observing that /3 2 is the variance of the 
operator (dHo/dp x ), which is always positive and van- 
ishes solely when the initial field distribution ip(x, y, 0) is 
an eigenfunction of (8Hq/ ' dp x ), i.e. of p x = —id x . Since 
such eigenfunctions are delocalized plane waves, it then 
follows that the variance of (dHo/dp x ) is strictly positive 
for any initial beam distribution carrying a finite power, 
regardless of the specific form of Hq. 



C. Self-collimation regime 

Beam self-collimation (also referred to as beam sub- 
diffraction) is a well known phenomenon occuring in 
homogeneous arrays and, more generally, in photonic 
crystal structures with engineered band structure i?o(p) 
showing points of local flatness (see, for instance, [20j ) . 
The simplest example of sub diffraction is the 'arrest' of 
beam spreading in a one-dimensional tight-binding lat- 
tice with nearest-neighboring couplings, which was ob- 
served in Ref.[|| using relatively broad Gaussian beams 
at an incidence angle set in correspondence of an inflex- 
ion point of the band dispersion curve. Though it is well 
understood that in such a regime diffraction is cancelled 
solely at low orders, it was perhaps overlooked the fact 
that self-collimation does not modify the beam broaden- 
ing scaling law [Eq.(17)]. In other words, self-collimation 
will correspond to a reduction of the coefficient f3 2 for 
special initial field distributions, but not to a change of 
the scaling law describing beam broadening. If we con- 
sider, for the sake of simplicity, a one-dimensional lat- 
tice and assume that the spectrum F(k) of the exciting 
beam, defined as F(k) = 1/(2tt) / dxi(j(x, 0) exp(— ikx), 
is narrow at around its mean fco, the value of /? 2 , as 
given by Eq.(19), can be expanded in series of moments 
Il=J dk(k-k y\F(k)\ 2 {I =2,3,4,...) as 



Pi = b 2 I 2 + b 2 b 3 h + ^(Ab 2 b i + ib 2 3 ) U- h ll 2 + ... (20) 



where bi is the value of the derivative (d l H /dk l ) evalu- 
ated at k = k . As rapidly goes to zero as the order I 
increases, Eq.(20) shows that at the points k of the dis- 
persion curve where 62 (and possibly 63, 64, ...) vanishes 
beam broadening is reduced. We will refer to such points, 
where the dispersion curve Ho(k) is locally flat, to as self- 
collimation points [note that the condition iJg(fco) = is 
not requested]. 

As an example, let us consider the simplest one- 
dimensional waveguide array in the nearest-neighboring 
approximation, considered in Ref. [3] to demonstrate self- 
collimation effects. The Hamiltonian Hq has the form 
Hq = — 2Acos(p), and the self-collimation points are lo- 
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cated at p — ±7r/2. From Eq.(19) one obtains 



Pi = 2A 2 



A 



1 - RC ^2 < C "+2^ 



^c* (c n+ i - c„_i) 



(21) 



For a bell-shaped (e.g. Gaussian-shaped) and flat beam 
incident onto the array at a given tilting angle 9 (nor- 
malized to the Bragg angle), we may write c n — 
\c n \ exp(— iirdn), and one obtains 

Pl{9) = 2A 2 [l - k\ + (k\ - k 2 ) cos(2tt0)] (22) 

where k\ and k 2 are defined by 



= ^2 i c « c «+ii > «2 = ^2 



CnC n +2\ 



(23) 



Generally, it turns out that k 2 > k 2 , so that the minimum 
of j3 x is attained at 9 — ±1/2, i.e. at the self-collimation 
points as expected from Eq.(20). Conversely, the max- 
imal diffraction (maximum value of /3 X ) is attained at 
normal incidence (9 = 0). The ratio between the mini- 
mum and maximum values of (3 X , given by 



= 1/2) 
Px(6 = 0) 



1 + K2 



2k\ 



1 — K2 



(24) 



may get very small for a broad input beam. To illustrate 
this point, let us consider as an example the following 
beam distribution at the input plane : |c„| = TVa'™', 
where the parameter a (0 < a < 1) determines the spot 
size of the input beam (a —> for single waveguide ex- 
citation, and a — > 1 for a plane wave excitation), and 
Af = [(1 — a 2 ) /(l + a 2 )] 1 / 2 is a normalization factor. For 
such a field distribution, the values of coefficients K\ and 
n 2 can be evaluated in a closed form, and read 



2a 



1 + a 2 



K2 



a 2 (3 



1 + a 2 



(25) 



The ratio T between the diffraction parameters at subd- 
iffractive (9 = 1/2) and normal incidence (9 = 0) regimes 
takes then the form [see Eq.(24)] T = [(1 - a 2 ) /(I + 
a 2 )] 1 / 2 . Note that, for a very broad beam excitation 
(a —> 1), both Ki and k 2 gets close to 1, f3 x tends to 
vanish [see Eq.(22)], and the diffraction length ~ l/j3 x 
diverges independently of beam tilting angle 9, as ex- 
pected for a very broad input beam. However, in this 
case the ratio of diffraction lengths in the normal (9 = 0) 
and subdiffractivc (9 = 1/2) regimes, which scales as 
~ r, tends to vanish as T ~ (1 — a) 1 / 2 . Conversely, for a 
very narrow input beam (a — > 0), both k\ and k 2 vanish 
and the diffraction length ~ 1/Ac turns out to be inde- 
pendent of tilting angle and given by ~ 1/(V2A) [see 
Eq.(22)] as expected for single waveguide excitation. 



D. Discrete far-field pattern and anomalous beam 
decay 

In spite of the fact that the asymptotic law describing 
beam broadening due to diffraction is the same for dis- 
crete and continuous media, a deep difference is found 
when analyzing the decay behavior of the field intensity 
versus propagation distance and the far-field diffraction 
patterns. For the sake of simplicity, we will consider the 
diffraction problem in one dimension, though the results 
may be generalized to the two-dimensional diffraction 
problem. We then rewrite Eq.(2) as 



id z tp{x, z) = H (p)ip(x, z), 



(26) 



where p = p x = —id x . For the usual paraxial one- 
dimensional diffraction problem in a homogeneous con- 
tinuous medium, one has Hg(p) = p 2 /2, whereas for dis- 
crete diffraction in a one-dimensional waveguide array 
one has H (—p) = H (p) (~tt < p < tt) and H' Q (p) = 
at p — and at the band edges p = ±tt. The most 
general solution to Eq.(26) can be written as 

tp(x, z) = J dkF(k) exp[ikx - iH (k)z] (27) 

where the spectrum F(k) is determined by the beam dis- 
tribution at the input plane ip(x, 0) 



F(k) = — I dxtp(x, 0) exp(— ikx) 
2tt J 



(28) 



(J dx — > J2n' x ~ * n an d = n ) c n in the discrete 
diffraction problem). Our aim is to calculate the decay 
behavior of the field amplitude ip( x J z ) as the propaga- 
tion distance z increases, at either a constant x position 
(for instance at x = 0) or along a line x = az, where a 
is a constant parameter defining the 'observation angle' 
of the diffracted pattern. Note that, as the observation 
angle a is varied, the function ipo(a;z) = ij)(x = az,z) 
corresponds, for large values of z, to the far field diffrac- 
tion pattern. We need thus to calculate the asymptotic 
behavior of the integral 



iftoia; z) — I dkF(k) exp[izg(k)} 



for 



where we have set 



g(k) = ak-H (k). 



(29) 



(30) 



For this purpose, we may use the methods of station- 
ary phase or steepest descend (see, for instance, (24j). 
which predict that the asymptotic behavior of ^o(a; z) 
as z — > oo depends on the existence and of the order of 
stationary points of the phase g(k) inside the integration 
domain. 

Let us first consider the continuous diffraction problem, 
Hq(p) = p 2 /2, and re-derive the well-known result that 
the amplitude iPq(oc;z) decays as ~ \^fz for any obser- 
vation angle a and the far-field pattern is proportional 
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FIG. 1: (color online) Beam propagation in a one-dimensional 
tight binding lattice with nearest neighboring coupling terms, 
showing far-field properties of discrete diffraction. In (a) 
the intensity distributions \ip(x, z)\ 2 are plotted, in arbitrary 
units, for propagation distances z = 0, z = 5/A, z — 10/A, 
z = 15/A, 2 = 20/A, and z = 25/A, where A is the coupling 
rate between adjacent waveguides. The inset in (a) shows the 
Gaussian spectrum F(k) of the beam (wk = 1.5). In (c) the 
intensity distribution \ip(x,z)\ 2 is depicted for a propagation 
distance z = 20/A as numerically calculated by Eq.(27) (solid 
curve) and by the approximate relation (33) (dotted curve, al- 
most overlapped with the solid one). In (c) the beam intensity 
is plotted at a propagation distance z = 500/ A, clearly show- 
ing the dominance of two peaks at the diffraction cone edges 
(self-collimation directions) and the onset of three different 
decay laws at \a>\ > 2A, a — ±2A, and \a\ < 2A. 



to the Fourier spectrum of the input (near-field) distri- 
bution. In this case, g(k) = ak — k 2 /2 has a unique 
saddle point at k = fco = ct, with g"(ko) = —1^0; 
therefore, provided that the spectrum F(k) has a nonva- 
nishing component at k = ko and F(ko) does not diverge 
[25} , according to the method of steepest descend one has 



ip {a; z) ~ F(a)d — exp[iza J /2 - in /4] (31) 

as z — > 00. From Eq.(31) we obtain the well-known 
result of paraxial diffraction theory that the amplitude 
tpo(a; z) of the beam decays as ~ Xyfz for any observation 
angle a (25[, and that the far- field diffraction pattern is 
shaped as the Fourier spectrum F(a) of the near-field 
distribution. This scaling law may be referred to as the 
normal scaling law, in the sense that the beam intensity 



I oc I -0| 2 decays as ~ X/z whereas the beam spot size w x 
increases asymptotically as ~ z [see Eq.(17)], the product 
Iw x being constant according to the power conservation 
law. 

For the discrete diffraction problem, we prove now that 
the decay law is generally slower than ~ X/yfz at the 
observation angles corresponding to self-collimation, and 
that the far-field pattern is peaked at such angles and 
does not reproduce the spectrum F of the near-field dis- 
tribution. To this aim, let us observe that, according 
to the steepest descend method, the slowest decay term 
in the integral of Eq.(29) comes from the saddle points 
3' (fco) = of largest order. In particular, if fc is a saddle 
point of order n > 2, i.e. g(k) ~ g(ko) + [g^ n \ko)/n\\(k — 
fco)' 1 for fc close to fco (g( n >(ko) 7^ 0), the contribution of 
the saddle point to the integral in Eq.(29) for large values 
of z is given by [24| 



F(k ) 



|z3 (n) (fco)| 1/f 



-(n! 



n 



exp 



izg(k )±i — 
2n 

(32) 

Therefore, the decay law for ipo(a; z) scales as ~ z 
where n is the highest order of the saddle points of g(k), 
provided that F(fco) 7^ 0. In the case of diffraction in 
a homogeneous continuous medium, the order of saddle 
point is always n — 2. To determine n for the discrete 
diffraction problem, let us note that the dispersion curve 
Ho(k) admits of at least a couple of inflection points, say 
at fc = ±fc , at which ffg(fco) = 0. These points cor- 
respond to the self-collimation directions introduced in 
Sec. II. C. Since g'(k) — a — H^k), the inflection points 
fc = ±fco turn out to be also saddle points when the ob- 
servation angle a is chosen equal to ff (±fco). Therefore, 
for the discrete diffraction problem the largest order n 
of saddle points is at least n = 3, and the decay law 
of ij)a{oL;z), at the two observation angles a — ±iJg(fco) 
corresponding to the self-collimation directions ±fc , is 
slowed down -as compared to continuous diffraction- to 
(at least) ~ z -1 / 3 . More generally, if the dispersion 
curve Ho(k) of the lattice is engineered to achieve a very 
flat behavior near a self-collimation point fc = fco, with 
g"(k ) = g'"(k ) = ... - g^Hh) = and g^(k Q ) ± 
(n > 3), the decay law of ipo(a;z) scales as ~ z 1 l n 
at the observation angle a = H'(ko). This scaling law 
of beam decaying in the discrete diffraction problem is 
therefore anomalous, in the sense that along the self- 
collimation directions the intensity decays slower that 
X/z, i.e. of the characteristic decay law that one might 
expect from power conservation arguments. This seem- 
ingly paradoxical circumstance may be solved by observ- 
ing that, for an observation angle a different than any 
of the self-collimation directions, the decay of t{jo(a;z) 
may be either normal (i.e. ~ Xj^fz) or even faster. More 
precisely, for a fixed value of a in modulus larger than 
ctmax = maxfc|.ff (fc)|, the function glk) given by Eq.(30) 
does not have saddle points on the real axis, and tpo(ce; z) 
decays as ~ X/z. Conversely, for \a\ < a max the equation 
g'(k) = a — H (k) =0 admits of at least one solution, 
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with g"(k) / for a second-order saddle point. In this 
case, according to the method of stationary phase the 
asymptotic behavior of ipo(a; z) for large values of z fol- 
lows the normal law ~ To summarize, tpo(a; z) 
scales: as ~ F(ko)z~ 1 / n at a self collimation direction 
a, where H (ko) — a and kg is a saddle point of order 
n > 3; as ~ F(ko)z~ 1 for an observation angle a outside 
the 'diffraction cone' \a\ > a max \ as ~ F(fc )z -1 / 2 in- 
side the diffraction cone region \a\ < a max but far from 
a self collimation direction. The far-field pattern of dis- 
crete diffraction tends therefore to confine light inside 
the diffraction cone \a\ < a max with asymptotic peaks 
at the propagation directions corresponding to the angles 
of self-collimation. 

This very general behavior may be illustrated more in 
details for the case of a tight-binding lattice in the 



J 



nearest-neighbor approximation considered in Sec. II. C, 
for which Ho(k) = — 2Acosfc. In this lattice model one 
has g'(k) = a — 2Asinfc, g"(k) — — 2Acosfc, so that the 
angle of diffraction cone is given by a max — 2A. Two 
saddle points of second-order are found at k = kg = ±ir/2 
for the observation angles a = ±a max , i.e at the edge of 
the diffraction cone, at which the far-field discrete diffrac- 
tion pattern is thus expected to show two peaks. For an 
observation angle a strictly inside the diffraction cone 
(\a\ < 2A), the equation g'{k) = has two solutions 
which are saddle points of first order since g"(k) ^ 0. 
The main contribution to the integral on the right hand 
side of Eq.(29) comes from these two saddle points, and 
can be calculated by the method of stationary phase, 
yielding explicitly 



ip(a;z) 



ip(a; z) 



Y ZsJA 2 - (a/2) 2 
(a > 0) 



V z^/A 2 - (a/2) 2 
(a < 0) 



{— iF(k(j) exp [iakQZ + 2iAz cos fco] + F( n ~ ko) exp [iaz(n — ko) — 2iAz cos fco]} 

{—iF(—k ) exp [— iak z + 2iAz cos fco] + F(—Tt + ko) exp [iaz(— n + fc ) — 2iAz cos fco]} 

(33) 



r 



where fco is the solution to the equation sin fco = |a/2A| 
in the interval < fco < 7r/2. It should be noted that 
the far-field discrete diffraction pattern given by Eq.(33) 
holds for \a\ < 2A. As \a\ approaches 2A from below, 
the two saddle points of second order coalesce into a sin- 
gle saddle point of third order, and this explain the di- 
vergence of Eq.(33) as \a\ — > 2A, i.e at the self colli- 
mation directions, where the decay is slower and scales 
as ~ z -1 / 3 . For |a| > 2A, i.e. outside the diffraction 
cone, there are not saddle points on the real axis and 
the decay is faster and scales as ~ 1/z. An example 
of far-field discrete diffraction for a beam with a Gaus- 
sian spectrum F(k) ~ exp[— (k/wk) 2 ] (— 7r < fc < n) 
is shown in Fig. I. In Fig. I (a) the intensity distribution 
\ip{x, z)\ 2 , as obtained by accurate numerical computa- 
tion of the integral entering in Eq.(27), is plotted for a 
few propagation distances z. For the sake of readabil- 
ity, at each propagation distance z the field intensity has 
been normalized to its peak value. The diffraction cone 
and the emergence of two intensity peaks at the self- 
collimation directions a = ±a max are clearly visible just 
after a propagation distance z of ~ 10 — 20 times the cou- 
pling length I/A [Fig. 1(a)]. Inside the diffraction cone, 
the intensity distribution at such propagation distances 
is very well fitted by the analytical far-field distribution 
given by Eq.(33), as shown in Fig. 1(b). At much longer 
propagation distances, the self-collimation peaks become 



dominant, and the appearance of three different scaling 
laws of beam decay (fast decay outside the diffraction 
cone |a;| > 2Az; normal decay inside the diffraction cone 
|x| < 2Az; slower decay at the self-collimation directions 
x = ±2Az) is very clearly visible, as shown in Fig. 1(c). 



III. BEAM PROPAGATION IN CURVED 
WAVEGUIDE ARRAYS 

Discrete diffraction of light waves in linear optical 
waveguide arrays can be controlled by introducing trans- 
verse index gradients or local phase slips, which may 
produce a kind of refocusing or re-imaging of beam dis- 
tributions along the propagation distances (see, for in- 
stance, [1, d, [IH [H, HH) similarly to what happens to 
light propagating in continuous lensguide media. In par- 
ticular, waveguide arrays with a curved axis provide a 
particularly interesting set up to manage discrete diffrac- 
tion for both monochromatic and polychromatic light 
SEMI- It is therefore of major interest to have gen- 
eral laws describing the global behavior of beam prop- 
agation in curved waveguide arrays. In addition, it is 
well known that for the problem of paraxial diffraction 
in homogeneous continuous media or, more generally, of 
paraxial propagation in elementary optical systems and 
lensguides, one can introduce special families of field dis- 
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tributions (such as the Gaussian beams) that propagate 
maintaining unchanged their functional shape (shape- 
invariant beams), and that field propagation may be sim- 
ply described by means of algebraic equations ruling out 
the evolution of some complex- valued beam parameters 
(such as the complex g-parameter for Gaussian beams; 
see, for instance, [23|)- A natural question is whether one 
can similarly introduce shape-invariant discrete beams, 
i.e. field distributions that do not change their functional 
shape when propagating along curved waveguide arrays. 
As the problem of discrete diffraction in waveguide arrays 
with curved axis or transverselly-imposed index gradi- 
ents is analogous to the problem of one-dimensional or 
two-dimensional Bloch oscillations of electrons in peri- 
odic potentials with an applied electric field or of cold 
atoms in optical lattices, some results are already avail- 
able in the literature. In particular, in recent works [28l — 
|30| an algebraic approach has been developed, capable of 
providing rather general results for wave packet center of 
mass evolution and wave packet spreading in certain lat- 
tice models. In this approach, after the introduction of a 
dynamical Lie algebra, an explicit form of the evolution 
operator is first derived, and then the expectation val- 
ues of operators are calculated in the Heisenberg picture. 
However, the question of existence of shape-invariant dis- 
crete beams and of their propagation in curved waveg- 
uide arrays does not seem to have been addressed yet. In 
this section, we present a generalization of Eqs.(12) and 
(17) describing the evolution of beam center of mass and 
beam width in curved waveguide arrays using the method 
of moments. Though similar results have been previ- 
ously published in Refs. [28l-[30j using an algebraic oper- 
ator approach, they are here re-derived for the sake of 
completeness using the method of moments, which does 
not require the explicit calculation of the evolution op- 
erator and the formulation of the problem in terms of a 
Lie algebra. In the subsequent section a family of shape- 
invariant discrete beams will be introduced, proving that 
their propagation in a generally-curved waveguide array 
is simply described by the evolution of a complex-g beam 
parameter, which plays an analogous role of e.g. the 
complex-^ parameter of Gaussian beams propagating in 
paraxial continuous optical systems. 
Let us consider monochromatic light propagation in a 
two-dimensional waveguide array with a curved axis de- 
scribed by the parametric equations x = xq(z) and 
y = yo(z); the coupled mode equations describing light 
transfer among coupled waveguides in the single-band 
and tight-binding approximations are an extension of 
Eq.(l) to include fictitious transverse index gradients in- 
duced by waveguide curvature and read explicitly 



verse index gradient entering in Eq.(34) may be also re- 
alized by applying a thermal gradient, or may describe 
lumped phase gradients (2(| or an abrupt tilt of waveg- 
uide axis direction Q, in which cases £ (z) shows a delta- 
like behavior. After introduction of a continuous function 
i/>(r, z) such that V(r„ >m , z) = c n ^ m (z), one can read- 
ily check that the discrete diffraction equations (34) are 
equivalent to the following continuous Hamiltonian prob- 
lem 



id z ip(r,z) = H(r,p)ip(r,z) 
(p = — iV r ) with Hamiltonian 

H = H (p)-£(z)-v, 



(35) 



(36) 



where Hq is the Hamiltonian of the homogeneous ar- 
ray defined by Eq.(3). The laws governing the evolution 
of beam center of mass and beam variance can be ob- 
tained by extending the method of moments described 
in Sec. II. B for the free diffraction problem. In general, 
the cascade of equations that one obtains by applying 
the Eherenfest equation (9) to (r), (x 2 ), (y 2 ) - and to the 
commutators found throughout the calculations- turns 
out to be unlimited for a general form of Ho, and a 
closed set of equations are found solely for special forms 
of Hq. Such a special circumstance is encountered in 
case of a one-dimensional waveguide array in the nearest- 
neighboring approximation, and in case of a rectangular- 
lattice waveguide array neglecting diagonal interactions. 
The first model corresponds to the Hamiltonian 



H(x,p) = — 2A cosp — £ x (z)x, 



(37) 



where A is the coupling rate between adjacent waveg- 
uides, and p — p x — —id x . The second model, which 
has been for instance considered in the experiment of 
Ref . [32j . is described by the Hamiltonian 

H(r,p) = ~2A X cos(p x )-2A y cas(p y ) -£ x (t)x-£ y (z)y, 

. (38) 

where A^ (A y ) is the coupling rate between adjacent 
horizontal (vertical) waveguides of the lattice [31 j - 



(34) 



A. One-dimensional array 

Application of the moment method to the one- 
dimensional Hamiltonian model (37) yields a set of closed 
coupled equations for the expectation values of operators 
x, 9, and of a; 2 , p and a, where 



where £(t) = £ x (t)u x + £ y (t)u y , £ x (z) = -(n s /X)x (z), 
£ y (z) — — (n s /X)y Q (z), n s is the refractive index of the 
waveguide substrate, and X — X/(2n) is the reduced 
wavelength of light. It should be noticed that the trans- 



6 = exp(ip) (39) 
P = \ {1 ~ exp(-2 ip )} (40) 
a = i {xexp(— ip) + exp(— ip)x} . (41) 
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Successive application of the Ehrenfest equation (9) 
yields the following equations for (x) and (9) 



d(x) 
dz 

dz 



= 2 Aim ((6»)) 



(42) 
(43) 



and the following coupled equations for (a; 2 ), (p) and (a) 
d(x 2 ) 



dz 
did 

dz 
dip) 

dz 



2ARe((<r)) (44) 
-2i£ x {z)(p)+i£ x {z) (45) 
4A(p)-i£ x (z)(a). (46) 



Equation (43) can be readily integrated, yielding the fol- 
lowing evolution equation for the beam center of mass 



(x(z)} = (ar(0))+2Im{<zofi*(z)} 
where we have set 



Sl(z) = / d£Aexp[-i<j>(£)], 
Jo 

(j>{ z ) = f de^(o, 

Jo 

q = ^2c* n (0)c n+1 (0). 

n 

Similarly, integration of Eqs.(45) and (46) yields 

(p(z)} = exp[— 2i(j){z)\ x 

x |(p(0)) + iexp[2^)]-i 

(cr(z)) = exp[— i<p{z)\ x 
1 



(47) 

(48) 

(49) 
(50) 



x { (<r(0)) + 4fi(») 



(p(0)) 



2fi*(z) 



(51) 



(52) 



Taking into account that A exp[— z0(z)] = dQ/dz and 
that 2Re{/*d£fi*(£)(dfi/d£)} = N(z)| 2 , substitution of 
Eq.(52) into Eq.(44) yields 

(x 2 (z)} = (x 2 (0)) + 2\Q{z)\ 2 + 2Re {qMz) - q 2 tf{z)} , 



where we have set 

q x = i^(2n-l)<(0)c„_ 1 (0) 

q 2 = 5>;(o) c „_ 2 (o). 

The beam size w x is then given by 



w 



.{z) = ^(x 2 (z)) - (x(z)y 



(53) 

(54) 
(55) 

(56) 



For a given field distribution c„(0) at the input plane, 
the evolution of the beam center of mass (x(z)) and beam 



size w x (z) are thus ruled by Eqs.(47), (53) and (56). Note 
that beam evolution depend on the input beam parame- 
ters qi, qi and qz -defined by Eqs.(50),(54) and (55)- and 
by the complex amplitude Q, (z), defined by Eqs. (48-49) 
and accounting for bending of waveguide axis. Note also 
that for straight arrays Q(z) = Az, and one thus retrieves 
the results of discrete diffraction derived in Sec. II. B, in 
particular the linear asymptotic increase of w x with z. 
The condition for diffraction suppression, i.e. a non- 
secular growth of w x (z) with z, is that Q(z) remains a 
limited function of z as z increases. This condition is 
always satisfied for a constant value of £ x , which corre- 
sponds to circularly-curved waveguides and to the onset 
of the optical analogue of Bloch oscillations [t| . Similarly, 
for periodic axis bending with spatial period A, £ x (z) is 
a periodic function of z, and the condition of boundness 
of Q(z) is given by 



/ ^exp[-#(0] = 0, 
Jo 



(57) 



which is precisely the condition of 'dynamic localization' 
previously investigated in Refs.fl^. fl3|. 



B. Two-dimensional array 

For the two-dimensional waveguide array model (38), 
the moment equations turn out to decouple into two set of 
equations, similar to Eqs. (42-46), separately acting onto 
the x and y directions. The evolution equations for the 
beam center of mass (x), (y) are then given by 



(x(z)) - 

(y(z)) - 

where we have set 

4>x,y{z) : 

and 



(x(0))+2lm{q Qx n x (z)} 
(y(0)) + 2Im{q 0y a* y (z)} 



d£A XtV exp[-i<j> Xt y(£)] 



d£,£x, y {Q 



<lo. 



qo v 



x = 2] C »,m( ) c n+l,m( () ) 



(0). 



(58) 
(59) 



(60) 
(61) 

(62) 
(63) 



Similarly, the beam sizes w x and w y , defined as 

w x (z) = ^(x 2 (z)) - (x(z)) 2 (64) 

(65) 



v y (z) = ^(y 2 (z)) - (y(z)) 2 , 



are calculated using Eqs. (58-59) and the following evolu- 
tion equations for (x 2 (z)) and (y 2 (z)) 

(x 2 (z)) = (x 2 (0))+2\n x \ 2 + 2Re{q lx n x -q 2x n 2 x } (66) 
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(y 2 (z)) = (y 2 (0)) + 2\n y \ 2 + 2Re {q ly Q y - q 2y n 2 y ) (67) 
where we have set 









q-ix 




(69) 


qi y 


= z^(2m-l)< m (0)c„, m _ 1 (0) 


(70) 


Q2y 




(71) 



IV. SHAPE-INVARIANT DISCRETE BEAMS 

The existence of shape-invariant beams, i.e. families of 
field distributions that propagate without changing their 
functional shape, is well-known for paraxial propagation 
in Gaussian optics or in continuous lensguide media (see, 
for instance, [23]). Here we address the related problem 
of investigating the existence of shape-invariant discrete 
beams, i.e. field distributions that do not change their 
functional shape when propagating along waveguide ar- 
rays with arbitrarily curved optical axis. This is a rather 
challenging problem because no general method capable 
of constructing shape-invariant beams seems to be avail- 
able. However, for the simple waveguide array models 
considered in the previous section, a family of shape- 
invariant discrete beams can be introduced in a simple 
manner. Owing to their functional form, such beams are 
referred to as discrete Bessel beams. 



A. Discrete Bessel beams in one-dimensional 
arrays 

Let us consider a one-dimensional waveguide array 
with an arbitrarily curved optical axis. In the tight- 
binding and nearest neighboring coupling approxima- 
tions, light propagation is described by the following set 
of coupled-mode equations 

ic n = -A(c n+ i + c„_i) - nf(z)c n (72) 

where f(z) describes the rate of transverse index gradi- 
ent induced by waveguide bending lumped waveg- 
uide tilting 3] or locally imposed phase changes among 
adjacent waveguides [26j as discussed previously. Let us 
fist observe that, if c n (z) is a solution to Eq.(72) corre- 
sponding to a given initial field distribution c n (0), then 
for an arbitrary integer no 

g n (z) = c„_„ (z)exp|m J j ( 73 ) 

is the solution to Eq.(72) corresponding to the trans- 
lated initial field distribution g n (0) = c n _„ (0). There- 



fore, apart from an unimportant phase change, shape- 
invariant beams remain invariant for an arbitrary trans- 
verse translation on the lattice. 

Let us tentatively search for a solution to Eq.(72) of the 
form 

Cn{z) — J n (a) exp(-icm) (74) 

where J n is the Bessel function of first kind of order n, 
and a = a(z), a = cr(z) are unknown functions which de- 
pend on propagation distance z, but not on lattice site n. 
Note that, as £„n|J„(a)| 2 = and [£ n \J n {a)\ 2 n 2 } = 
a 2 /2, the parameter a is related to the beam size w x 
[Eq.(6)] by the simple relation w x — a/v2, whereas a 
defines a transverse tilt of the beam 'phase front'. Sub- 
stitution of Eq.(74) into Eq.(72) and taking into account 
the identities of Bessel functions J„+i(a) + J„_i(a;) = 
(2n/a)J n (a) and J„_i(a) — J n +i(a) = 2J! n (a), one ob- 
tains that Eq.(74) is indeed a solution to Eq.(72) pro- 
vided that a and a satisfy the coupled equations 

a = -2Asincr (75) 
2A 

a = coser— /. (76) 

a 

Owing to the functional form of c„, we will refer such 
shape-invariant beams to as discrete Bessel beams. Let 
us define a complex-g parameter for the discrete Bessel 
beam (74) according to 

q{z) = a(z) exp[ia(z)] (77) 

so that the modulus of the complex q parameter gives 
the beam spot size at propagation distance z, whereas 
its phase corresponds to the phase front gradient. From 
Eqs.(75) and (76) one readily obtains for the complex q 
parameter the following simple evolution equation 

^=-2iA-if(z)q. (78) 

The general solution to Eq.(78), for a given initial value 
q(0) at the z — input plane, is given by 

q(z) = e X p[-i<P(z)]{q(0)-2i^ d£A exp j (79) 

where 

m= f duio- (so) 

Jo 

The propagation of a discrete Bessel beam along a 
curved waveguide array is thus reduced to the prop- 
agation of its complex q parameter, which plays an 
analogous role of the complex-q parameter for Gaussian 
beams in lensguide media. The propagation law of the q 
parameter admits of a simple geometrical interpretation 
in the complex q plane. According to Eq.(78), for an 
infinitesimal propagation distance 5z the change of q{z) 
is given by the superposition of the two paths AB and 



(a) (b) 




FIG. 2: (a) Geometric construction of the evolution of the 
complex q parameter for an infinitesimal propagation distance 
Sz. The length of the segment AB is 2A8z, whereas the ro- 
tation angle is £7 = f(z)Sz. The points A and C correspond 
to q(z) and q(z + Sz), respectively, (b) Geometric representa- 
tion of a self-imaging array: the path followed by the complex 
parameter q(z), starting from the origin O, is closed. 
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Re( 9 ) 



FIG. 3: (color online) (a) Schematic of a one-dimensional 
waveguide array with a tilt of waveguide axis at z = zo, and 
the transformation induced on the complex q parameter by 
the tilt (inset), (b) Principle of beam collimation via waveg- 
uide axis tilt [tilt angle 9 — X/(4an s )], and path followed by 
the complex q parameter for single waveguide input excitation 
from z = to z = d + (inset). 



BC shown in Fig.2(a). The path AB, of length 2SzA, 
accounts for discrete diffraction and corresponds to a 
change of q(z) along the imaginary q axis; the path BC 
is due to the transverse index gradient which produces a 
clockwise rotation by the angle Sj — f(z)6z around the 
origin O of the complex plane. It is interesting to note 
that, since J„(0) = 8 n< o, for q(0) = the discrete Bessel 
beam (74) reduces to the well-known impulse response 
of a tight-binding array with nearest neighbor couplings 
(see, for instance, [33(). 

To appreciate the usefulness of the ^-parameter descrip- 
tion and some properties of discrete Bessel beams, let us 
now discuss a few examples and applications. 

Propagation of discrete Bessel beams in homoge- 
neous arrays. 

For a homogeneous array (/ = 0), the propagation law 
of the complex-g parameter is simply given by 



q(z) = q(0) - 2iAz. 



(81) 



If we assume, for the sake of definiteness, that at the 
input plane z = the phase front of the beam is fiat, i.e. 
q(0) = a(0) = cvo real valued, the following propagation 
laws for beam size a and beam phase tilt a are derived 



(82) 
(83) 



From Eq.(82) we may introduce, as for Gaussian beams 
propagating in free space [2l} , the Rayleigh range zr and 
divergence angle 9<i such that at(zR) — \/2a>o and 9d = 
lim 2 ^oo a(z)/z, i.e. 




d 



2A 
2A. 



(84) 
(85) 



It should be noted that, as opposed to the case of 
Gaussian beams in free space -for which the Rayleigh 
range zr is proportional to the square of the spot size 
ao at the beam waist and the diffraction angle 9^ is 
inversely proportional to ao- for discrete Bessel beams 
the Rayleigh range zr is proportional to the spot size 
ao at the beam waist whereas the divergence angle is 
independent of the beam spot size and always equal to 
the diffraction cone angle introduced in Sec. II. D. This 
peculiar property is closely related to the very general 
result, proven in Sec. II. D, that the far field of discrete 
diffraction in a homogenous waveguide array is peaked 
at the observation angles corresponding to the flattest 
points (self-collimation points) of the band dispersion 



Transformation of a discrete Bessel beam through a 
waveguide axis tilt. A tilt of the waveguide axis at z = zq 
by a (small) angle 9 corresponds to impressing a phase 
shift 
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2tt 



■a9n* 



(86) 



between adjacent waveguides, where a is the waveguide 
spacing and n s the effective index of propagating modes 
[see Fig. 3(a)]. Light propagation across the tilt can be 
thus modelled by assuming f(z) = jS(z — z ) in Eq.(72), 
and its effect on the complex q parameter is to produce 
a rotation around the origin of the complex plane by an 
angle 7 [see the inset of Fig. 3(a)]. 

A tilt of the waveguide axis may be used to 'collimate' 
a discrete beam, as schematically shown in Fig. 3(b). 
Here a single waveguide is initially excited at the input 
plane, and after a propagation distance d the axis of 
the array is tilted by an angle 9 = \/{Aan s ) such that 
7 = 7r/2. The 90° rotation of the q parameter in the 
complex plane due to axis bending [see the inset in 
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Fig. 3(b)] brings the q parameter on the real axis, with 
a zero phase gradient a = and an enlarged beam size 
a = 2Ad. The axis tilt thus plays a similar role of a 
collimating lens for a diverging Gaussian beam. Note 
however that, contrary to a conventional lens, the tilting 
angle 9 to achieve beam collimation is independent of 
the distance d between source point (at z = 0) and 
the lens plane (z = d). Figure 4 shows an example 
of beam collimation in a 6-cm-long one-dimensional 
array as obtained by numerical analysis of the scalar 
wave equation for the electric field envelope E(x, z) 
propagating in the structure based on a standard beam 
propagation method. Figure 4(a) shows a pseudocolor 
map of the intensity beam evolution \E(x, z)\ 2 along 
the structure when a single waveguide is excited in its 
fundamental mode at the input plane z = and the 
waveguide axis is tilted at a distance d — 2 cm from the 
input plane [horizontal dotted curve in Fig. 4(a)]. The 
refractive index profile n(x) used in the simulations is 
depicted in Fig. 4(b), and the values of other parameters 
are A = 1.55 /im, n s — 1.52, and a = 11 /im, corre- 
sponding to a tilting angle 9 — \/(Aan s ) ~ 23.2 mrad. 
For the sake of readability, the intensity distribution 
is plotted with the waveguide axis z unfolded along a 
straight line. Note that the numerical results provide 
a realistic behavior of beam propagation beyond the 
couple-mode equation approximation, accounting for 
radiation losses and coupling to higher-order bands due 
to axis bending. These latter effects, however, are very 
small for the parameter values adopted in the simu- 
lations, and the coupled-mode equation model works fine. 

A geometric interpretation of the self-imaging condition 
and polygonal Block oscillations. An array of length 
d shows a self-imaging property (also referred to as 
diffraction cancellation or dynamic localization) , when- 
ever |c n (<i)| 2 = |c„(0)| 2 for any initial field distribution. 
The dynamic localization condition has a rather simple 
geometric interpretation in the complex q plane. In fact, 
if the array is excited in waveguide n = 0, q(0) = and 
to achieve self-imaging after a propagation distance d 
one has necessarily to have q(d) = q(0) = 0, i.e the path 
described by the complex q parameter, starting from 
the origin O of the complex plane, should be closed 
[see Fig. 2(b)]. Owing to the translational invariance 
of discrete Bessel beams [Eq.(73)], this condition is 
also sufficient. From Eq.(79), the closed-path condition 
q(d) = q{0) = yields 



dz exp[i<j>(z)] = 



(87) 



which is precisely the condition for dynamic localization 
derived originally by Dunlap and Kenkre in Ref. [331] . 
An application of the geometric condition of dynamic lo- 
calization is that of polygonal Bloch oscillations. Let us 
consider a waveguide array whose axis forms an (open) 
polygonal curve of large (mean) radius R made of a se- 
quence of straight segments of same length 6 and with tilt 
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FIG. 4: (color online) (a) Pseudocolor map of beam intensity 
propagation in a waveguide array with one axis tilting as ob- 
tained by numerical simulations, showing beam collimation. 
(b) Refractive index profile of the waveguide array used in 
the numerical simulations. The values of other parameters 
are given in the text. 



angle 9, as shown in Fig. 5(a). The function <fi(z), defined 
by Eq.(80), is thus a staircase function, which increases in 
steps of 7 = (27T /X)a9n s [see Eq.(86)] at z = 6,26,36,... 
(the coordinate z is measured along the polygonal curve). 
After a propagation d = (N + 1)6 from the input z = 
plane, where N is an integer number, it then follows that 

/ dz exp[i<j>(z)] = 6^^ exp(?7n). 

J n=0 

The sum of complex numbers (phasors) on the right hand 
side of Eq. (88) can be done analytically and has a well- 
known geometric interpretation; in particular, if 7 satis- 
fies the condition 7 = 2tt/ (N + 1), i.e. if the tilt angle 9 
is given by 



A 



an s (N+l) 



(89) 



the sum on the right hand side of Eq.(88) vanishes, and 
the condition for self-imaging is attained. An example of 
the self-imaging property of a polygonal waveguide array 
is shown in Fig. 5(b) for the case N = 5. The figure 
depicts a characteristic breathing mode corresponding to 
a single waveguide excitation at the input plane. The 
waveguide array parameters are the same as in Fig. 4, 
and a sequence of axis tilts are placed at distances 6=1 
cm one to the next. The tilt angle 9, chosen according to 
Eq.(89), is 9 ~ 15.5 mrad, yielding a self-imaging plane 
at d = (N + 1)6 = 6 cm, as clearly shown in Fig. 5(b). 
Note that, in the limit 6 — > 0, N — > 00 and b/9 -> R 
finite, the polygonal of Fig. 5(a) approximates an arc of 
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FIG. 5: (color online) (a) Schematic of a polygonal waveg- 
uide array for the observation of Bloch oscillations, (b) Pseu- 
docolor image of beam intensity propagation in a 8-cm-long 
polygonal array showing a Bloch oscillation breathing mode. 
The refractive index profile of the waveguide array used in the 
numerical simulations is the same as in Fig. 4(b). The values 
of other parameters are given in the text. 

a circumference of radius R, and the condition (89) for 
self-imaging is satisfied for a propagation distance 

d=(N + l)b^— (90) 
n s a 

which is the spatial period of Bloch oscillations on a 
curved waveguide array (radius of curvature R) previ- 
ously considered in Refs.|9l. Hoj|. The usual Bloch oscil- 
lations on a curved waveguide array may be therefore 
viewed as a limiting case of Bloch oscillations on a polyg- 
onal array. 

B. Discrete Bessel beams in two-dimensional arrays 

A simple extension of the analysis of Sec. IV. A can be 
done for a two-dimensional rectangular-lattice waveguide 
array with nearest-neighboring coupling when the diago- 
nal coupling is neglected. This model is described by the 



coupled mode equations 

iCn,m — ^x(Cn+l,m Cn — l,m) ^y(Cn,m+l H~ Cn,m — l) 

- f x (z)ncn im - f y (z)mcn,m (91) 

where f x , y (z) describe the rates of transverse index gra- 
dients induced by waveguide bending or lumped waveg- 
uide axis tilting along the x and y directions. Since Eqs. 
(91) admit of separable solutions c„, m (z) = c n (z)c m (z), 
with c n (z) and c m (z) solutions to the one-dimensional 
problem (72) with A = A. x ,y and f(z) — f x ,y(z), a two- 
dimensional discrete Bessel beam has the form 

c n , m (z) = J n {a x ) Jm(ot y ) exp(-ier x n - a y m). (92) 

The complex-g parameters of the beam along the x and 
y directions are defined by 

q x (z) = a x (z) exp[i(x x (z)] , q y (z) = a y {z) exp[ia y {z)} 

(93) 

and their evolution is ruled out by the equations 

Qx,y ^i^oc.y ^f X ,y\Z) (94) 

which have a similar geometric interpretation as that dis- 
cussed in Sec. IV. A. The propagation properties of two- 
dimensional discrete Bessel beams in homogeneous ar- 
rays, across tilted axis regions or polygonal curves are the 
same as those investigated for one-dimensional beams, 
and arc therefore not further discussed here. 



V. CONCLUSIONS 

In this work, a comprehensive study of discrete diffrac- 
tion and linear propagation of light in homogeneous and 
curved waveguide arrays has been presented. In partic- 
ular, general laws describing beam spreading, beam de- 
cay and discrete far-field patterns in homogeneous arrays 
have been derived using the method of moments and the 
steepest descend method, and some remarks on the well- 
known self-collimation regime have been pointed out. In 
curved arrays and within the nearest neighboring cou- 
pling approximation, the method of moments has been 
extended to describe the evolution of global beam pa- 
rameters. This method provides an alternative means to 
algebraic operator techniques recently proposed in other 
physical contexts to study general properties of Bloch 
oscillations (28l - [30| . Finally, a family of shape- invariant 
discrete beams -referred to as discrete Bessel beams ow- 
ing to their functional form- has been introduced. It has 
been shown that propagation of such beams in curved 
waveguide arrays is simply described by the evolution of 
a complex q parameter, which plays a similar role to the 
complex q parameter used for Gaussian beams in contin- 
uous lensguide media. A few applications of the q param- 
eter formalism are discussed, including beam collimation 
via waveguide axis tilting, a geometric interpretation of 
the self-imaging effect in waveguide arrays, and optical 
Bloch oscillations on a polygonal array. 
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